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Preface 


The  purpose  of  this  study  was  to  investigate  the 
feasibility  of  using  the  two-parameter  negative 
exponential  distribution,  with  scale' set  equal  to  one,  as 
a  baseline  for  a  robust  estimator  for  the  location 
parameter  (minimum  life)  of  selected  life  distributions. 
This  estimator  could  be  used  in  determining  cost  and 
performance  estimates  for  systems  containing  many 
electronic  components  (e.g.,  satellites),  where  expense 
and  small  quantities  make  it  impractical  to  conduct 
enough  tests  to  determine  underlying  time  to  failure 
distributions  with  certainty. 

A  computer  model  was  built  using  Monte  Carlo 
techniques  to  generate  time  to  failure  data  for  several 
underlying  distributions.  Five  estimates  for  the  minimum 
life  were  computed  from  this  data  and  compared  to  the 
actual  minimum  life.  The  best  estimator  of  the  five  being 
compared  was  shown  to  be  the  minimum  variance,  unbiased, 
maximum  likelihood  estimator  for  the  minimum  life  from 
Che  two-parameter  negative  exponential  distribution. 
However,  the  limits  of  the  usefulness  of  this  estimator 
still  need  to  be  determined. 

The  development  of  this  thesis  required  a  great  deal 
of  help  from  others.  I  would  like  to  thank  my  thesis 
advisor.  Dr.  Albert  H.  Moore,  for  his  guidance  and 


technical  assistance.  My  limited  theoretical  knowledge 
could  have  made  this  thesis  an  impossible  task.  His 
patient  explanations  of  underlying  theory  and  its 
relationship  to  the  real  world  quite  often  lifted  a  very 
dense  fog.  I  would  also  like  to  thank  my  reader.  Dr. 
Joseph  T.  Cain.  His  feedback  concerning  the  written 
portion  of  my  thesis  was  invaluable. 

Finally,  I  would  like  to  thank  my  fellow  classmates, 
the  finest  group  of  people  I  have  ever  worked  with.  Their 
support  and  camaraderie  over  the  last  18  months  helped 
ease  the  way  through  a  difficult  program.  In  particular  I 
would  like  to  thank:  Major  Dennis  Charek  and  Captains 
Karen  Barland,  Jim  Porter,  and  John  Sours  for  their  help 
in  finding  background  information  and  debugging  my 
computer  program;  Captains  Jim  Martin  and  Norm  Jarvis  for 
their  moral  support  at  my  thesis  defense  -  the  morning 
after  Thanksgiving;  and  First  Lieutenant  Stacy  Brodzik 
for  her  support  and  encouragement,  and  for  letting  me 
know  when  I  was  taking  things  too  seriously. 


Linda  M.  Allen 


Table  of  Contents 


Page 


Preface  . 11 

List  of  Figures  .  vl 

List  of  Tables  . vll 

Abstract  .  vlll 

I.  Introduction  . 1 

Problem  .  2 

Objectives  .  2 

Scope  .  2 

II.  Background  . 4 

Life  Distributions  .  4 

The  Exponential  Distribution  .  6 

The  Welbull  Distribution  .  7 

The  Gamma  Distribution  .  8 

The  Log-normal  Distribution  . .  9 

Robust  Estimation  . 11 

Minimum  Distance  Estimation  .  12 

The  Kolmogorov  Statistics  . .  12 

The  Cramer-von  Mises  Statistic  .  13 

The  Anderson-Darling  Statistic  .  13 

Measures  of  Effectiveness  .  13 

Mean  Square  Error  . 13 

Relative  Efficiency  . 14 

III.  Methodology  . 15 

Overview  of  Method  . . , .  15 

Random  Number  Generation  .  15 

Parameter  Estimation  . 16 

Maximum  Likelihood  Estimation  .  16 

Distance  Estimation  .  16 

Investigation  of  Robustness  of  Estimator  ....  16 

IV.  Results  .  18 

V.  Conclusions  and  Recommendations  .  20 

Appendix  A;  Graphs  of  Representative  Probability 

Density  Functions  . 21 

Appendix  B:  Tables  I-IV,  Mean  Square  Errors  for 

Location  Parameter  Estimates  .  26 


i  V 


List  of  Tables 


Table 

I.  Mean  Square  Errors  for  Location  Parameter 

Estimates  (Sample  size=8)  . 

II.  Mean  Square  Errors  for  Location  Parameter 

Estimates  (Sample  size=12)  . 

III.  Mean  Square  Errors  for  Location  Parameter 

Estimates  (Sample  size  =  16)  . . . 

IV.  Mean  Square  Errors  for  Location  Parameter 

Estimates  (Sample  size=20)  . 

V.  Relative  Efficiencies  of  Location  Parameter 

Estimates  (Sample  size=8)  . 

VI.  Relative  Efficiencies  of  Location  Parameter 

Estimates  (Sample  size=12)  . 

VII.  Relative  Efficiencies  of  Location  Parameter 

Estimates  (Sample  size=16)  . . . 

VIII.  Relative  Efficiencies  of  Location  Parameter 

Estimates  (Sample  size=20)  . 


This  investigation  determined  that  for  time  to 
failure  distributions  that  are  moderate  deviations  from 
the  negative  exponential  distribution,  a  robust  estimate 
of  the  minimum  life  could  be  arrived  at  by  assuming  the 
underlying  distribution  was  exponential  and  using  the 
minimum  variance,  unbiased,  maximum  likelihood  estimator. 
It  was  found  that  estimators  using  the  minimum  distance 
statistics  of  Kolmogorov,  Cramer-von  Mises,  and 
Anderson-Darling  did  not  perform  well  with  the  asymmetric 
distributions  explored  in  this  thesis.  However,  they  may 
still  prove  useful  for  life  distributions  with  larger 
shape  parameters  (i.e.,  for  distributions  that  are  not 
"moderate"  deviations  from  the  negative  exponential). 

The  analysis  was  accomplished  by  using  Monte  Carlo 
techniques  to  generate  random  samples  of  time  to  failure 
data  from  specific  distributions,  and  using  this 
empirical  data  to  estimate  the  actual  minimum  life  of  the 
distribution.  Five  estimators  were  explored:  the  minimum 
variance,  unbiased,  maximum  likelihood  estimator  of  the 
two-parameter  negative  exponential  distribution;  the 
first  ordered  statistic;  and  the  three  minimum  distance 
methods  (the  theoretical  distribution  was  assumed 
negative  exponential).  The  performance  of  these 
estimators  was  evaluated  by  comparing  their  mean  square 


Robust  Estimation  of  the  Location  Parameter  of  Life 

Distributions 

I .  Introduction 

Space  assets  have  been  playing  an  ever  greater  role 
in  the  defense  of  this  nation.  Satellites  can  be  placed 
at  altitudes  that  give  them  a  world  view.  Unfortunately, 
once  in  place  these  satellites  are  mostly  inaccessible  to 
us  -  if  they  break,  they  are  expensive  to  repair,  if  they 
can  be  reached  to  repair  at  all.  In  addition,  the  space 
environment  is  hostile,  which  tends  to  increase  the 
failure  rate, 

A  typical  satellite  consists  of  hundreds  of  parts, 
all  of  which  have  to  be  made  to  exacting  standards.  There 
is  a  great  deal  of  expense  involved  in  the  designing, 
building,  testing,  and  placement  of  these  satellites.  In 
order  to  make  decisions  regarding  the  number  of 
satellites  to  be  built,  or  whether  they  should  be  built 
at  all,  there  must  be  some  way  of  estimating  the  expected 
minimum  life  of  the  satellite. 

Testing  different  satellite  components  can  produce  a 
rough  time  to  failure  data  base  that  can  be  used  to 
estimate  minimum  life.  In  the  case  of  a  satellite,  a  very 
expensive  piece  of  equipment,  it  is  hoped  that  the 
expected  minimum  life  is  greater  than  zero. 


Satellites  are  generally  expensive  and  manufactured 
in  small  quantities.  It  may  be  difficult  to  collect 
enough  time  to  failure  data  to  identify  the  underlying 
distribution  with  certainty.  What  is  desirable  is  an 
estimator  which  gives  a  useable  value  in  most  cases  even 
when  the  assumption  about  the  underlying  distribution  is 
incorrect.  Such  an  estimator  is  termed  "robust"  (Parr  and 
Schucany,  1979). 

PROBLEM 

A  robust  estimation  method  is  desired  for  the 
location  parameter  of  life  distributions.  This  location 
parameter  is  the  guaranteed  minimum  life. 

OBJECTIVES 

The  overall  objective  of  this  thesis  was  to  develop 
a  robust  estimator  for  the  location  parameter  of  life 
distributions,  using  the  two-parameter  negative 
exponential  distribution,  with  scale  equal  to  one,  as  the 
baseline , 

SCOPE 

Only  four  families  of  distributions  were  considered: 
exponential,  Weibull,  gamma,  and  log-normal.  The  sample 
sizes  from  each  distribution  were  small.  Four  sample 
sizes  were  used:  8,  12,  16,  and  20. 


Three  minimum  distance  methods  were  investigated: 
Kolmogorov,  Cramer-von  Mises,  and  Anderson-Darling.  The 
robustness  of  the  estimators  will  be  evaluated  using  two 
measures  of  effectiveness:  mean  square  error  and 
relative  efficiency. 


I I .  Background 


If  the  underlying  distribution  of  a  specific 
population  is  known,  predictions  can  be  made  regarding 
samples  from  the  population.  For  example,  if  it  is  known 
that  the  time  between  failures  of  the  components  of  a 
particular  type  of  satellite  is  exponentially  distributed 
then  predictions  can  be  made  concerning  the  performance 
of  that  satellite  over  a  specified  period  of  time. 

Unfortunately,  the  underlying  distribution  of  a 
population  is  usually  something  that  has  to  be  determined 
from  a  relatively  small  sample  from  that  distribution. 
There  is  a  great  deal  of  literature  concerning  the 
determination  of  underlying  distributions  and  the 
statistical  inferences  that  may  be  made  from  these 
determinations.  Miller's  and  Freund's  text  on  probability 
and  statistics,  and  Kapur's  and  Lamberson's  text  on 
reliability  engineering  were  used  extensively  in  the 
development  of  this  thesis. 

Life  Distributions 

This  thesis  is  primarily  concerned  with  life 
distributions.  When  electronic  components  are  tested  two 
of  the  more  important  statistics  collected  are  time  to 
failure  and  time  between  failures. 

Time  to  failure  has  been  modeled  with  several 
distributions,  the  exponential,  gamma,  Weibull,  and 


log-normal  being  of  particular  interest  to  this  thesis 
(Banks  and  Carson,  1984:134).  The  exponential 
distribution  is  most  commonly  used  for  life  testing 
applications,  with  the  Weibull  distribution  probably  the 
second  most  common  (Kapur  and  Lamberson,  1977:233,  291). 

Banks  and  Carson  had  the  following  to  say  about  time 
to  failure: 

"If  only  random  failures  occur,  the  time  to 
failure  distribution  may  be  modeled  as 
exponential.  The  gamma  distribution  arises  from 
modeling  standby  redundancy  where  each 
component  has  an  exponential  time  to  failure... 

When  there  are  a  number  of  components  in  a 
system  and  failure  is  due  to  the  most  serious 
of  a  large  number  of  defects,  or  possible 
defects,  the  Weibull  distribution  seems  to  do 
particularly  well  as  a  model...  The  log-normal 
distribution  has  been  found  to  be  applicable  in 
describing  time  to  failure  for  some  types  of 
components,  and  the  literature  seems  to 
indicate  increased  use  of  this  distribution  in 
reliability  models."  (Banks  and  Carson, 

1984:134) 

The  gamma  and  Weibull  distributions  both  include  the 
exponential  distribution  as  a  special  case  (Banks  and 
Carson,  1984:132). 

The  exponential  distribution  is  often  used  in  life 
testing  applications  because  it  is  easy  to  apply  (Kapur 
and  Lamberson,  1977:233).  Also,  in  the  case  of  a  total 
system  composed  of  many  components  with  different  failure 
distributions,  the  time  to  system  failure  distribution 
will  approach  the  exponential  (Kapur  and  Lamberson, 
1977:236).  Of  particular  interest  in  this  thesis  is  the 
two-parameter  exponential  distribution.  This  distribution 


form  is  used  in  nonzero  minimum  life  situations,  i.e., 


situations  in  which  there  is  an  initial  period  of  no 
failures  (Kapur  and  Lamberson,  1977:258). 

Following  is  additional  information  concerning  the 
four  distributions  of  interest. 

The  Exponential  Distribution.  The  two-parameter 
exponential  distribution  has  a  probability  density 
function  given  by 


f (x;m,d) 


(1/m)  exp{ -(  x-d ) /ml  ,  xi  d>0,  m>0 


0  ,  elsewhere 


( 


where  m  equals  the  mean,  and  the  parameter  d  is  the 
minimum  life.  The  mean  of  this  probability  density 
function  is  m+d  (Kapur  and  Lamberson,  1977:258). 

The  cumulative  distribution  function  is  then  given 
by 


F(x) 


'o 

(1/m)  exp{ - ( t-d ) /m}  dt 

0 

=  I  -  exp{ -( x-d ) /m) 


,  X  <  d 

( 

,  x2  d  >0 ,  m> 


In  the  situation  where  n  items  are  placed  on  life 
test,  and  the  test  is  terminated  at  the  time  of  the  rth 
failure,  minimum  variance,  unbiased,  maximum  likelihood 


estimators  for  m  and  d  are  m  and  d  respectively  and  are 
defined  by 


m'  =  ^  (n-r)(x^  -x^) 

(T=T) 


(3) 


and 


d'=x^-(m'/n)  (4) 

where  x^  is  the  first  ordered  statistic  (Kapur  and 
Lamberson,  1977:258-259).  This  thesis  primarily 
considered  the  case  where  n*r,  i.e.,  uncensored  samples. 

The  Weibull  Distribution.  The  cumulative 
distribution  function  for  the  three-parameter  Weibull 
distribution  is  given  as 


' 

F(x;b,a,d)  = 


0  ,  X  <d 

(5) 

l-exp{-((x-d)/(b-d) )®}  ,  x>d 


where  a>0,  b>0,  and  d20.  a  is  the  shape  parameter  or  the 
Weibull  slope,  b  is  the  scale  parameter  or  the 
characteristic  life,  and  d  is  the  location  parameter  or 
the  minimum  life  (Kapur  and  Lamberson,  1977:292). 


The  Gamma  Distribution.  The  probability  density 


function  for  the  gamma  distribution  is 


f  (  X  ;  a  ,  b  ) 


{  1  /  (  b®  r(  a )  )}  (  X  ^  ^)exp  {- (  x-d  )  /  b  } 

for  x>  0  ,  a>0,  b>  0 

(6 


0 


,  elsewhere 


where  a  is  the  shape  parameter,  b  is  the  scale  parameter 
and  r(a)  .is  a  value  of  the  gamma  function  defined  by 


r  (a) 


a-1  -X  , 
x  e  d  X 


(7 


Through  integration  by  parts,  the  above  equation 
reduces  to 


r(a)  =  (a-l)r(a-l)  (8 

for  a>0,  and 

r(a)  =  r(a-l)!  (9 

when  a  is  a  positive  integer.  (Miller  and  Freund, 
1977:117)  The  exponential  distribution  is  a  special  casi 
of  the  gamma  distribution  when  a  is  equal  to  one. 


The  cumulative  distribution  function  is  given  by 


F(t)  = 


,  t<0 


(10 


{l/(b^r(a))}  t^  ^  exp{ -(t-d)/b}  dt , 


OS  tSx 


If  a  is  an  integer,  successive  integration  by  parts 
yields 


F(x)  =  I^^g{(l/b)x}  exp{ -(x-d)/b} 


(11) 


k! 


(Kapur  and  Lamberson,  1977:25). 

The  Log-normal  Distribution.  "The  log-normal 
distribution  occurs  in  practice  whenever  we  encounter  a 
random  variable  which  is  such  that  its  logarithm  has  a 
normal  distribution."  (Miller  and  Freund,  1977:114) 

The  probability  density  function  for  the  log-normal 
distribution  is  given  by 

{1/(0  t/Ii )  } 

exp{-(l/2)((ln  t  -u  )/a)  },  t>0 


where  -oo<u<oo  and  O>0  .  Here,  the  logarithm  of  the  random 
variable  t  has  a  normal  distribution,  i.e.,  a  random 
variable  x  defined  as  x=ln  t  will  be  normally  distributed 
with  a  mean  of  U  and  a  standard  deviation  of  0(Kapur  and 
Lamberson,  1977:19).  The  mean  of  the  log-normal 
distr  bution  is 

2 

E(t)  =  exp{u  +  0  /2}  (13) 

and  the  variance  is 

■  2  2 

V(t)  =  }  {e"^  -  1}  (14) 


(Kapur  and  Lamberson,  1977:20). 

The  log-normal  distribution  has  a  cumulative 
distributon  function  of 


0 


T<0 


F(t) 


(15) 


'^(  l/(Ta/2i)) 

°  exp{-(l/2)((ln  T  -u)/0^)}dT,  OSTSt 


(Kapur  and  Lamberson,  1977:20). 


Robust  Estimation 


The  parameters  of  a  suspected  underlying 
distribution  can  be  estimated  using  the  information 
provided  by  a  random  sample  from  that  population. 
Statistical  methods  that  produce  estimates  that  are 
relatively  insensitive  to  assumptions  about  the 
underlying  distribution  have  been  termed  robust  methods 
(Crow  and  Siddiqui,  1967).  These  robust  parameter 
estimates  should  continue  to  perform  well  under  moderate 
deviations  from  the  suspected  distribution  (Parr  and 
Schucany,  1979:2). 

In  this  thesis  the  two-parameter  negative 
exponential  distribution  was  used  as  a  baseline,  with  the 
Weibull,  gamma,  and  log-normal  distributions  providing 
the  "moderate  deviations".  These  distributions  were 
chosen  because  of  their  similarity  to  each  other  (see 
Appendix  A).  Their  probability  density  functions  all 
contain  an  exponential  factor,  and,  for  small  shape 
parameters,  they  have  similar  curves  when  graphed.  With 
small  sample  sizes,  like  those  used  in  this  thesis,  it 
may  be  possible  to  fit  the  data  with  all  four 
distributions,  and  impossible  to  say  with  certainty  which 
distribution  is  the  actual  underlying  distribution. 


Minimum  Distance  Estimation 

Minimum  distance  estimation  has  been  considered  as  a 
method  for  deriving  robust  estimators  (Parr  and  Schucany, 
1979  and  1982).  The  "distance"  referred  to  here  is  a 
discrepancy  measure  between  an  empirical  distribution 
function  and  a  theoretical  cumulative  distribution 
function  (Parr  and  Schucany,  1979:3). 

The  theoretical  cumulative  distribution  function 
used  in  this  thesis  is  given  in  equation  (2).  The  mean, 
m,  will  be  estimated  using  equation  (3),  and  the  minimum 
life,  d,  will  be  initially  estimated  using  equation  (4). 
The  F(x)  values  so  computed  are  the  values  used  in  the 
distance  statistics.  The  initial  value  for  the  minimum 
life  is  then  varied  through  several  iterations  of  the 
distance  statistics  until  the  distance  is  minimized.  The 
value  of  the  minimum  life  that  achieves  this  minimum 
distance  is  recorded  as  an  estimate  for  the  actual 
m i n imum  life. 

Three  distance  statistics  were  tried  in  this  thesis: 
the  Kolmogorov,  Cramer-von  Mises,  and  Anderson-Darling. 
The  computing  formulas  for  these  statistics  are  given 
below;  their  theoretical  development  is  beyond  the  scope 
of  this  thesis. 

The  Kolmogorov  Statistics. 

UiSn 


D  =  max 


(17) 


D  =  max  ( ,  D“  ) 

The  Cramer-von  Mises  Statistic . 


(18) 


W"=  =  (2i-l)/2n)^  +  (l/12n)  (19) 

The  Anderson-Darling  Statistic. 

A"=  -  (2i-l)(ln  2.+  ln(l-Zj^^^_.  ))}  /n  -  n  (20) 

(Stephens,  1974:730-731). 

Measures  of  Effectiveness 


Mean  Square  Error.  The  mean  square  error  was 
computed  using  the  following  formula: 

MSE  =  I®^i(cl  -  d')"^  (21) 

s 

where  d  was  the  actual  minimum  life  of  the  underlying 
distributions  the  samples  were  generated  from,  d'  was  the 
estimated  minimum  life,  and  s  was  the  number  of  samples 
generated.  The  sample  sizes  were  varied  (8,  12,  16,  and 
20),  and  the  number  of  samples  taken  for  each 
configuration  of  input  parameters  and  sample  size 
remained  constant  at  one  thousand. 

A  small  mean  square  error  would  indicate  that  the 
estimated  minimum  life  deviated  little  from  the  actual 
minimum  life,  while  a  large  mean  square  error  would 
indicate  that  it  deviated  to  a  greater  extent. 


computed  as  follows: 


Relative  =  MSE  of  "best”  estimator  (22) 

Efficiency  MSE  of  comparison  estimator 

For  the  underlying  exponential  distribution,  the 
best  estimator  for  the  minimum  life  was  chosen  to  be  d’ 
(see  equation  (4)).  For  the  other  three  distributions, 
the  first  ordered  statistic  (x^)  was  chosen  as  the  best 
estimator,  primarily  because  it  was  easy  to  obtain.  A 
relative  efficiency  greater  than  one  indicated  the 
comparison  estimator  performed  better  than  the  "best" 
estimator.  Five  different  minimum  life  estimates  were 
compared  in  this  thesis:  the  exponential  estimate  (d'; 
see  equation  (4)),  the  first  ordered  statistic,  and  the 
three  minimum  distance  estimates. 
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III.  Methodolo 


Overview  of  Method 

The  two-parameter  negative  exponential  distribution 
was  chosen  as  the  baseline  for  computing  the  minimum  life 
estimates.  Random  samples  were  taken  from  the 
exponential,  Weibull,  gamma,  and  log-normal 
distributions.  Sample  sizes  of  8,  12,  16,  and  20  were 
used.  Each  sample  was  assumed  to  be  from  an  exponential 
distribution,  and  the  mean  and  location  parameters  were 
estimated  using  equations  (3)  and  (4).  These  estimated 
parameters  were  then  used  in  the  theoretical  cumulative 
distribution  function  used  in  computing  the  Kolmogorov, 
Cramer-von  Mises,  and  Anderson-Darling  distances  (see 
equations  (2),  (18),  (19),  and  (20)). 

One  thousand  runs  were  made  for  each  configuration 
of  sample  size  and  input  parameters.  The  mean  square 
errors  and  relative  efficiencies  for  each  of  the  minimum 
life  estimates  were  then  computed. 

Random  Number  Generation 

The  Monte  Carlo  method  is  a  technique  that  uses 
random  or  pseudorandom  numbers  for  solution  of  a  model. 
There  are  arithmetic  codes  available  at  many  computer 
centers  for  generating  sequences  of  pseudorandom  digits, 
where  each  digit  (0  through  9)  occurs  with  approximately 


equal  likelihood  (Rubinstein,  1981:11). 


The  pseudorandom  number  generators  used  in  this 
thesis  were  routines  from  the  International  Mathematics 
and  Statistical  Libraries,  Inc.  (IMSL).  Four  routines 
were  used:  GGEXN  to  produce  exponential  random  deviates, 
GGWIB  to  produce  Weibull  random  deviates,  GGAMR  to 
produce  gamma  random  deviates,  and  GGNLG  to  produce 
log-normal  random  deviates.  Ten  was  added  to  each  of  the 
random  deviates  to  simulate  a  minimum  life  of  ten.  This 
value  was  chosen  arbitrarily. 

Parameter  Estimation 

Maximum  Likelihood  Estimation.  For  the  nonzero 
minimum  life  situation  the  mean  and  minimum  life  were 
estimated  using  equations  (3)  and  (4)  respectively.  The 
underlying  distribution  was  assumed  to  be  exponential. 

Distance  Estimation.  The  theoretical  cumulative 
distribution  function  was  assumed  to  be  exponential,  and 
the  maximum  likelihood  estimates  for  the  mean  and  minimum 
life  were  used  in  equation  (2)  in  the  computation  of  the 
values  used  in  the  distance  estimation  equations 
(equations  (16),  (17),  (19),  and  (20)). 

Investigation  of  Robustness  of  Estimator 

One  thousand  samples  from  each  configuration  of 
input  parameters  and  sample  size  were  taken,  and  five 


estimates  for  the  minimum  life  were  made.  The  mean, 
standard  deviation,  mean  square  error,  and  relative 
efficiency  were  then  computed.  The  results  are  presented 
in  table  form  (see  Appendices  B  and  C)  and  are  analyzed 
in  Chapter  4. 


IV .  Results 


Appendix  A  contains  Figures  1  through  4.  These  are 
representative  graphs  of  the  probability  density 
functions  of  the  four  families  of  distributions 
considered  in  this  thesis. 

Appendices  B  and  C  contain  tables  summarizing  the 
mean  square  errors  and  relative  efficiencies  of  the 
minimum  life  estimators.  There  were  five  estimators 
considered:  the  exponential  estimator  (d';  see  equation 

(4)),  the  first  ordered  statistic,  and  the  three  minimum 
distance  estimators  (Kolmogorov,  Cramer-von  Mises,  and 
Anderson-Darling)  . 

An  inspection  of  Tables  I  through  VIII  shows  that  in 
most  cases  the  exponential  estimator  (d')  was  superior, 
with  the  Anderson-Darling  estimator  being  next  best. 
However,  the  minimum  distance  estimators  require  much 
more  time  computationally  than  the  exponential  estimator. 
In  order  for  one  of  them  to  be  considered  for  robust 
estimation  purposes  it  would  have  to  perform  much  better 
than  the  exponential  estimator. 

The  first  ordered  statistic  was  superior  when  the 
shape  parameter  equaled  .5.  More  runs  with  shape 
parameters  less  than  one  would  be  needed  to  determine 
just  when  the  first  ordered  statistic  becomes  a  better 
estimator  than  the  exponential  estimator. 
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From  a  visual  inspection  of  the  graphs  in  Appendix  A 
it  may  be  argued  that  as  shape  parameters  become  greater 
than  one,  the  underlying  distribution  may  no  longer  be  a 
"moderate"  deviation  from  the  exponential  distribution. 

In  these  cases,  a  robust  estimator  based  on  the 
exponential  distribution  could  not  be  expected  to  produce 
good  results.  Many  more  runs  with  different  shape 
parameter  values  would  be  needed  to  better  define  what 
could  be  considered  a  moderate  deviation  from  the 
exponential  distribution. 

The  IMSL  routine  ZXMIN  was  used  to  minimize  the 
minimum  distance  functions.  It  was  noted  that  the 
Kolmogorov  function  was  often  terminated  due  to  rounding 
errors  before  it  reached  a  minimum  value.  The  Kolmogorov 
values  were  left  in  the  tables  for  comparison  purposes 
should  someone  develop  an  algorithm  that  would  truely 
minimize  the  function. 
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V .  Conclusions  and  Recommendations 


The  exponential  estimator  for  the  minimum  life  ( d '  ) 
may  provide  a  robust  estimate  for  the  minimum  life  of 
distributions  that  deviate  moderately  from  the 
exponential.  What  constitutes  a  moderate  deviation  still 
needs  to  be  determined.  Variations  in  the  shape  parameter 
greatly  affected  the  mean  square  errors  of  the 
estimators.  More  tests  would  have  to  be  made  in  order  to 
determine  the  boundaries  of  acceptable  shape  parameter 
values  for  this  particular  robust  estimation  method. 

The  use  of  histograms  might  help  with  determining  if 
the  shape  parameter  is  too  large  to  use  the  exponential 
estimator.  However,  the  small  sample  sizes  typically  used 
with  robust  estimation  may  not  provide  enough  data. 

The  minimum  distance  estimators  did  not  do  well  with 
the  asymmetric  distributions  used  in  this  tliesis.  Further 
study  using  a  modified  minimum  distance  method  might  be 
feasible,  but  tremendous  improvement  in  the  mean  square 
errors  would  have  to  be  achieved  in  order  to  justify  the 
much  greater  expense  for  the  computer  time  required. 


Figure  2  . 


PDFs  for  several  gamma  distributions  when  0=1 
(Banks  and  Carson,  1984:145) 


MEAN  SQUARE  ERRORS  FOR  LOCATION  PARAMETER  ESTIMATES 


MHAN  SQIIARK  HKUOKS  FOR  I.OCATION  FARAMKTER  ESTIMATES 


MEAN  SQUARE  ERRORS  FOR  LOCATION  PARAMETER  ESTIMATES 


( Mu= 1 ,  Si gma=2 ) 
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TABLES  V  -  V 


RKl.ATIVR  EFFICIENCIES  OF  LOCATION  PARAMETER  ESTIMATES 


(Mu= 1 ,  Si gma  =  2  ) 
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PROSRAM  GKU«2(GKUA2C«TAPei0=CKUA20> 

UNSERLrING  DISTRIBUTION  IS  GAMMA. 
UNCSNSORCO  SAMPLES. 

FOUR  SAMPLE  SIZES  -  8tl2«lS«  AND  20. 

SCALEsl.O 

SHAPE  =  2.a  AND  4.9 


PROGRAM  function: 

-  GENERATES  GAMMA  RANDOM  NUMBERS. 

-  SORTS  RANDOM  NUMBERS. 

-  SAWES  THE  FIRST  ORDERED  STATISTIC. 

-  ESTIMATES  SCALE  AND  LOCATION  PARAMETERS 

ASSUMING  exponential  DISTRIBUTION. 

-  REESTIMATES  LOCATION  PARAMETER  USING 

KOLMORGOROV  DISTANCE. 

-  REESTIMATES  LOCATION  PARAMETER  USING 

CRAMER-VON  MISES  DISTANCE. 

-  REESTIMATES  LOCATION  PARAMETER  USING 

ANOERSON-OARLING  DISTANCE. 

-  COMPUTES  MEANS.  STANDARD  DEVIATIONS.  MEAN 

SQUARE  ERRORS.  AND  RELATIVE  EFFICIENCIES. 


VARIABLES  USED: 

N  -  SAMPLE  SIZE. 

I.L.P'.O  -  LOOP  COUNTERS  AfiO  ARRAT  INDICES. 

NRUNS  -  NUMBER  OF  RUNS  tNUMBER  OF  SAMPLES >* 

I  OPT  -  ZXMIN  INPUT  -  OPTIONS  SELECTION. 

MAXFW  -  ZXMIM  INPUT  -  MAXIMUM  NUMBER  OF  FUNCTION 
EVALUATIONS  ALLOWED. 

NSIO  -  ZXMIN  INPUT  -  NUMBER  OF  DIGITS  ACCURACY 
REQUIRED. 

NPARAM  -  ZXMIN  INPUT  -  NUMBER  OF  PARAVETERS. 
lER  -  ZXMIN  OUTPUT  -  ERROR  PARAMETER. 

0  -  ACTUAL  MINIMUM  LIFE  (LOCATION  PARAMETER*. 

M  -  SCALE  PARAMETER. 

TEMP  -  TEMPORARY  VARIABLE. 

R<N*  -  RANOCM  NUM9ER  ARRAY. 

Z<NI  -  Z  VALUE  ARRAY  (CDF*. 

X(NPARAM»  -  ZXMIN  INPUT/OUTPUT  -  PARAMETER 
VALUE. 

XISAVE  -  PREVIOUS  X<NPARA.M>. 

F  -  ZXMIN  OUTPUT  -  VALUE  OF  FUNCTION  BEING 
MINIMIZED. 

F3AVE  -  PREVIOUS  F. 

Hdl.OdI  -  ZXMIN  OUTPUT  -  WORKING  VECTORS. 

W<3>  -  ZXMIN  OUTPUT  -  ZXMIN  OUTPUT  - 
VECTOR. 

U(l)  -  NORM  CF  THE  GRADIENT. 

W<2I  -  NUMBER  OF  FUNCTION  EVALUATIONS 

performed. 

■J<3)  -  AM  ESTIMATE  OF  THE  NUMBER  OF 
SIGNIFICANT  DIGITS  IN  THE  FINAL 
PARAMETER  ESTIMATE. 

SHAPE  -  SHAPE  PARAMETER  (USED  IN  GAMMA  ANC 
VEIBULL  OISTRIBUTIONS  ONLY*. 

UK(4''*  -  WORKING  VECTOR  (USED  IN  GAMMA 

distribution  only* . 

MU  -  MEAN  (USED  IN  LO'J-NERMAL  DISTR  IH.UT  ICN 
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I  !• 

SIOfiA  -  STANDAOD  QEtflATIOM  CUSED  IN  LOG-NCRMAL 

ozsTRiaurioN  only). 

£D<NRUNS>  -  INITIAL  ESTIMATE  CF  THE  LOCATION 
PARAHETERt  ASSUMING  AN  EXPONENTIAL 
OtSTRieUTION. 

EMINRUNSI  -  ESTIMATE  OF  THE  SCALE  PARAMETER# 
ASSUMING  AN  EXPONENTIAL  DISTRIBUTION. 

XKNRUNSl  -  FIRST  ORDERED  STATISTIC. 

EOKINRUNSI  -  KOLMOGOROV  LOCATION  PARAMETER 
ESTIMATE. 

EOVINRUNSt  •  CRAMER-VON  MISES  LOCATION 
PARAMETER  ESTIMATE. 

EOA(NRUNS>  -  ANDERSON-DARLING  LOCATICN 
PARAMETER  ESTIMATE. 

OPLUSiDMINUS  -  KOLMOGC»OV  MAXIMUM  POSITIVE 
AND  NEGATIVE  DEVIATIONS. 

OSUM.KSUMtVSUM,ASUM«MSUM,XSUM  -  INTERMEDIATE 
VALUES. 

XOtXK.XU.XAiXHtXX  -  SAMPLE  MEANS. 

SOtSK tSU.SA .SMtSX  -  SAMPLE  STANDARD 
CEVIATIONS. 

MSEO.MSEKtMSEUtMSEAtMSEMtMSEX  -  MEAN  SQUARE 
ERRORS. 

OSEEO  -  INPUT  SEED  TO  RANDCM  NUMBER 
GENERATORS. 


SUBROUTINES  USEOt 

SORT  -  aOBBLE  SORTS  RANDOM  NUMBERS#  PLACING 
smallest  number  in  RCl). 

GOEXN  -  (PROGRAM  EKWA  ONLYI  IMSL  ROUTINE  TO 
GENERATE  EXPONENTIAL  RANDCM  NUMBERS. 

GGUIB  -  (PROGRAMS  UKUAl  AND  UKUA2  ONLY)  IMSL 
ROUTINE  TO  GENERATE  VEIBULL  RANDOM 
NUMBERS. 

3GAMR  -  (PROGRAMS  GKVAl  AND  6KVA2  ONLY)  IMSL 
ROUTINE  TO  GENERATE  GAMMA  RANDOM 
NUMBERS. 

GGNLG  -  (PROGRAMS  LKHAl  AND  LXVA2  ONLY)  IMSL 
ROUTINE  TO  GENERATE  LOG-NORMAL 
RANDOM  NUMBERS. 


INPUT  files;  none 


OUTPUT  files: 

EKHAO  -  PROGRAM  EKUA 
GKMA’O  -  PROGRAM  CKUAl 
GKVAZO  -  PROGRAM  GKWAO 
WKUAIO  -  PROGRAM  WKUAl 
UKUA20  -  PROGRAM  WKMAO 
LXVAIO  -  PROGRAM  LMUAl 
LKMA20  -  PROGRAM  LKUA2 


COM  (ON  R|-g,P#':T,EM#»l#dDK#:Dl##:CA,x;SAVE.F3AVE 
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lOPTtHAXFNtMSICtNPARAMiIER 
<<EAL  OfntTE'«P,R(2C>  •2(2C>tXlSAUE«FSAUEt 
X(  tt  <K(1  )  tGCl  >  tFtMf3>t 
ton  JQC)  (EHt  !<)  001 1  Kill  COO)  • 
EDK(ICOC)  •EDUClCCClfCOAdOOC  >« 

□  PLUS  «0M  INUS  tSHAPE  tHU>SIGHA  tV)K  (AO  1 1 

□  SUHtKSUH<'ilSU»<ASUMtnSUH.  XSUMt 
X3«XK«  XWtXAt  XN*X«t 
SOiSKtSUtSAtSMtSXt 
MSEOfMEEK.RSEU .MSEA# HSEM, HSEX 

OOURLE  PRECISION  OSEED 

•*  UERSET  SUPPRESSES  TERHINAL  ERRORS. 

CALL  UERSEKO  tLEVOLD) 

DSEE0=12IAE7. CO? 

NRUNS  =  1300 

1  =  1.? 

0=13.0 

NPARAF=1 

*<SIG=A 

MAXFN=500 

IOPT=0 


**  TRO  SHAPES  (2.0  ANO  A.O) 


SHAPE=2.0 
DO  a  L=lf2 


•  •  RUN  A  SAHPLE  sizes 


N=8 

00  A  0=ltA 

•  •  MAKE  lOOC  RUNS 

00  n  P  =  ltNP,UNS 


••  GENERATE  GAMMA  RANDOM  NUMBERS 


CALL  GGA»«R<OSEED.SHAPE*.N.RK.R* 


••  AOO  OFFSET  FOR  MINIMUM  LIFE 


00  S  I=lfN 
R(  I  )  =  R<  ll»0 

continue 


••  SORT  RANDOM  NUMBERS 


CALL  SORTINtRI 


(JUOUU  UUUUUU  UUU<JOU  uuuouu 


X1(PI=R(1) 


**  ESTIMftTE  SC«LE  AND  LOCATION  PARAMETERS 


*su«=o.; 

00  S3  I=2<N 

XSUf1=XSUM*R«  I  J-RIl  I 
so  CONTINUE 

EN{PI=XSUM/FL0AT<  N-1> 
E0<P»=Ra  )-£M<P>/FLOATtK) 


••  REESTIMATE  MINIMUM  LIFE  USING  KOLMOGOROV 
OISTANCE. 


X<1>=E0<P  » 

CALL  2XMIN<  KOISTtNPARAMtNSICtMAXFNt  lOPTt 
X»H»G.F»U.1EP.> 

EOK<P>=X( 1) 


•  •  REESTIMATE  MINIMUM  LIFE  USING  CRAMER-VON 
MISES  OISTANCE. 


X<U=EO<PI 

XISAVEsRIl) 

CALL  ZXMIM  UOIST.NPARAMtNSIC.MAXFNtlOPT, 
X.H.G.f »W,IER> 


EOU(P>sX<l| 


*•  REESTIMATE  MINIMUM  LIFE  USING 
ANOERSON-OARLING  OISTAACE. 


xa>=EO(P> 

XiSAVE=R<l) 

FSAVE=.3RR 

CALL  ZXMr.(AOIST*NPARAP.t,SIG.MAXFN,ICPT, 
♦  Xf  Hf  G»  F f  f  lER  ) 

£0A(PI=X<1) 

10  CONTINUE 
C 
C 

C  ••  COMPUTE  SAMPLE  MEAN 

C 

C 

OSUM=0.0 
KSUM-0 .0 
MSU.N=1.0 
XSUM=;.0 

WSUM=J.I 
ASUM=0.3 
00  R8  1=1. ‘.RUNS 
0SUM  =  3S'JM*E0I  I» 

KSUM=KSUM.EOK  (I  I 

msum=msum»emi  I> 

XSU.'t  =  XSU.''-*Xl<I> 
mSU"  =  43UM«£D  .(<I  ) 

A3UM  =  A3'jri»EDA  <I  I 
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•r .  v«,  mT*  .  *. 


O'i  L  'J  «  1  I  .^uc 

'I0  =  3SU’* /FLOAT  tMFUNS  I 
XK  =  KSU'1/FLOAT(NI<UNS  » 

xh=msuh/float<nruns  » 

XX  =  XSUM/FLC*T  THRUMS  » 
XU=4SUH/FL0AT<riRUNSI 
XA  =  ASUH/FLOAT<NRUNS  > 


I  t* 


**  COMPUTE  sample  STANOARO  DEVIATIONS 


0SUM=0.0 
KSUM=0 .0 
MSUM=S.') 

XSUMsO.O 
USUM=0.0 
ASUMsO.O 
00  S9  I=l«NRUNS 

DSUM=OSUM*<EO<I >-XD> ••2.0 
KSUM  =  KSUM^(EOKt  I»-XK»^«2.I) 
.«SUM=M3UM*(  E.M  (1 1  -XM)  •*2  .0 
XSUM  =  XSUM*<  »1  <H-XXI^^2.0 
WSUM=WSUM*TEDU« I> -Xa>^^2.0 
ASUM  =  A3UfT»<EDAT  I»  -XA)^«2,3 
89  CONTINUE 

S0=<0SUM/FLCAT<NRUNS-1I  •••.5 
SKsTKSUM/FLOATTNRUMS-IT  /••.S 
SM=<MSUM/FLCAT«NR UNS-1)  )»».5 
SX  =  <XSUM/FL0AT<NRUNS-1>  »»».s 
SU=<USUM/FLOAT(MRUNS-l»  )«».5 
SA  =  <ASUM/FLOAr<MRONS-l  >  *»«,5 


•  •  COMPUTE  MEAN  SQUARE  ERRORS 


OSUMsC.'l 
XSUM=0.1 
MSUMsJ .0 
XSUMsO.O 
WSUM=0.: 

AsuM=o.': 

00  11  IrltNRUNS 

OSUM=OSUf*»<  0-E0«  I  I  •••2 
KSUH=KSUM»<0-E3K< :>!•• 
••SUM  =  MSUM*<  M-£M«  I  »»^»2 
XSUM=XSUM»<0-X1 «I l)^»2 
V3UM:v(SU"»<0-E3V<  !*»•■< 
ASUM  =  ASUM*<  C-EO  A<  I)  ••« 
11  CONTINUE 

MSEO=OSUM/FLOAT (NFUNS» 
MSEK=K3JM /FLOAT TMRUNS) 
MSEM=MSUM/[-LOAT(NRUNS» 
*iSEX=XSUM/FLOAT(NRU  NS> 
'1SEW=U3UM /FLOAT  (NRUNSI 
MSEA=A3UM/FL0AT«NRUNS) 


PRI'IT  TAHLE 


WRITEdO  t8C0) 

ROO  F0R''»T(  n** 'OUTPUT  FROM  GXWA2'//» 
VRITE(13t«;il 

801  FT9  matt  1  <  ,  "JNOfOLTI  NO  OISTRTSUTION  tS  OAMMA') 


•  -(i  I  -  I*.  »  =  .  .#  vru  .- 

805  FC^'^AT  <  IX  I  •LOCATIO'I  tO»  =«,F1C.5/ 

♦  1*. •scale  =*tFn.3/ 

♦  IXt  'SHAPE  =•  fFlC.3/ 

♦  Itf’SAMPLE  SIZE  <NI  ='tI3/ 

♦  IXt'NBR  OF  RUNS  (NRUNS>  =*»I5///) 

URITEdOiSlCI 

813  FORMATUBXt 'WEAN*  * 

URITEClO.aiS) 

8  13  FORNATClTVt  'SOUARE*  t9*  . 'RELATI  NE  •  .2  8*  i 

♦  •STA^OARO•) 

URITE(iat320> 

823  FORMAT < 18 »« 'ERROR 't  Tx,* efficiency 'flSX* 

♦  '•1EAN'«8X,'3EWIArT0N'l 

WRI TE 1 13 ««2S)MSEO»HS£X/MSEO*XO.SD 
325  FOR  iUTl  3X  I  'EXP*  t4t5X,F12.«l » 

XR ITE( 13 1830 IMSEK  tNSEX/nSEK t XK i SK 
8  3C  FORMAT*  3X.  'K0L'.«<5X.F12.«;  » 

wa  ITE«1C  .fl35IMSE8.HSEX/liSEytXX.SW 
335  FORMAT*  3X, 'CVH' ,A  (5XtFl2.8l  > 

yR  ITE<13 ,8NCIMSEA,MSEX/MSCA,XA,SA 
833  FORMAT*  4X  «•  AO' .**  5X  •F12.8>  > 
yR ITE* 13  tSASIMSEX tXXf SX 
8  33  F0RMAT*3X  t  'Xl'  t5X  tF  12.8«17Xt2*  5X,F12.8»  > 
yRI  rE*10  tSSOIMSEM.XH.SM 

853  FORMAT*  IXt 'SCALE'  t5XtF12.8tl7X,2*5**F12.8)) 
N=N*3 

8  CONTINUE 

SHAPE=SHAPE»2.0 
8  CONTINUE 
STOP 
END 


3URR0UTINE  S0RT(NtR> 

INTEGER  I.JtN 
REAL  rEMPiR<N) 

00  335  I=lfN-l 
00  310  Js!l*l.N 

IF*R<  JJ  .LT.  R<  I)  >  THEN 
C  ••  SUAP 

TEMP-:R*  I) 

R*I  l=<><  J) 

R*  JI^TEHP 
END  IF 
31C  CONTINUE 
335  CONTINUE 
RETURN 
ENO 


SUBROUTINE  KOI  ST*  NP  ARAM  ,x  tF  > 

CO “MON  RiNtPtEOtEMf  Xl *EDKtEDU» ECA«X IS A VE»FS AVE 
INTEGER  NPARAH.NtP 

REAL  X<NPARAM|  ,OPLUStOMI.NUS  .F,  TEMP,  XI SA  I/E,  FSA  VE  , 
»  Z*20 ) ,R*2C I ,E0*i:C3) ,EM*10CC ),CDK(1C00) , 

*  xmooci  ,EDU*1COC>,£OA*130C> 

C  ••  COMPUTE  Z*II 
00  41  1  =  1  , N 

Z<I )  =  :.C-EXP*-«R*  I)-X<NPAHAM) >/E“IP ) > 

SI  CONTINUE 
C  ••  FINO  MAX  0* 

OPLUS=3.0 
0  3  70  I=1,N 

TEMP  :FL3AT«  I)/FLOAT*N»-Z«  I» 

IF*TEMP  .OT.  OPLUSI  OFLUS  =  TEMP 
73  CONTINUE 
C  ••  FI  NO  max  0- 
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uni  uu  i-u  •  -j 
DO  T1  1=1, N 

TEMP  =  Z(  I)-FL0ATCI-1»/FL0ATCN» 

IFIT^IP  .GT.  OMINUS)  0P1NUS=TE:MP 
71  CONTINUE 

F=M4X<QPLUS,D MINUS! 

RETURN 

END 

SUBROUTINE  UDI S  T(  NP  AR  AP,X.F  ) 

COMMON  R,N,P,EOtEM,Xl,EDK«EOU,EOA,XlSAVE,FSAVE 
INTEGER  NPARAMtN.P 

REAL  X<NPARAM>,XSUH,F,TEMP,2<20),Rf20 ) , XIS A WE, FSA VE , 

*  XKIOOO)  ,EO(1000>  (EHdOOO  >  , 

*  EOK(igOO>,EOU<1000)tCOAflOOO> 

••  XlSAWE  IS  THE  PREVIOUS  VALUE  OF  X(NPARAM). 

ZXMIN  OCCASIONALLY  PUTS  AN  EXTREMELT  LARGE 
number  in  XTNPARAM).  mhen  this  happens, 

THE  PREVIOUS  VALUE  IS  USED. 

IF(X(NPARAM>  .LE.  R(l>>  THEN 
X1SAVE=X(NPARAM> 

ELSE 

X<NPARAM)=X1SAWE 

ENOIF 

•  •  COMPUTE  Z<I> 

DO  62  1  =  1  ,N 

2<I) =1. 0-EXP (-<R< I)-XtNPARAM) ) /EH(P) ) 

62  CONTINUE 

••  COMPUTE  CRAMER-VOM  MISES  DISTANCE 

XSUM=0 .0 
DO  SO  1  =  1, N 

TEMP=FL0AT<2*I-1I /FLCAT<2»N) 

XSUM  =  XSUM=x3UM»fZ  <I>-TEMP)«»2.0 
80  CONTINUE 

F  =  XSUM*FL0AT<1) /FLOAT  (12«N» 

RETURN 

END 

SUBROUTINE  AO  IST<  NPARAM,X  ,F  » 

COMMON  R,M,P,ED,EM,X1,E0K,ED«,EDA,X1SAVE,FSAVE 
INTEGER  NPARAM,N,P,FLAG 

real  x<nparam),f,z(20 j,r (20 >,xsu.m, temp, fsave, XI save, 

*  XKIOC:)  ,ED(1000)  ,EH(100C  I  , 

*  EOK< ICCC) ,EOU( 1000>,EDA( 1300) 

MINIMUM  life  not  greater  THAN  FIRST 
OROEREO  STATISTIC. 

IF{»(NPARAM)  .GT.  R«l)»  THEN 
»<NPARA  M)=X1SAVE 
ELSE 

x:SAVE=  <(NPARAM) 

ENOIF 

•  »  COMPUTE  Z<I» 

FLAG=0 

DO  63  1  =  1  ,N 

Z<n=l.:-EXP<-<R(  I)-X(NPARAMI  >  /EM(P  I  ) 

1F(  (ZC)  .LE.  C.)  .:R.  (Z(I>  .GE.  1.)  )FLAG  =  1 

63  CONTINUE 

IF  (FLAG  .£0.  ?)  THEN 

••  COMPUTE  ANCERSCN-CARLING  CISTANCE 
XS'JMrO.O 
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